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A COMPUTER PROGRAM FOR MAPPING SATELLITE-BORNE 


NARROW-BEAM ANTENNA FOOTPRINTS ON EARTH 


I. INTRODUCTION 

A computer program has been developed that computes the locus of 
intersection of a quadric cone and a sphere. The outputs of the program 
are a list of the longitude and latitude coordinates of the locus of 
intersection and a plot of the locus. It was written primarily to define 
the area of the earth covered by a narrow-beam antenna carried on a 
synchronous satellite in circular, near equatorial orbits. 

The program is basically an implementation of a report by S. L. 
Zolnay^ with some modifications added. The main modifications are the 
incorporation of an elliptical cross section antenna beam and computation 
of the beam vertex angles and corresponding locus of intersection for 
any signal level between 0.1 and 10 db down from the beam center. Any 
number of signals up to 10 may be input for each data set. The output 
will be plotted on the same graph. 

The program was written for use with a Cal Comp 570/563 off-line 
plotting system and uses the standard Cal Comp subroutines supplied with 
the system. The plot is drawn using linear longitude and latitude scales 
and non-linear scales such as Mercator scales cannot be used. 

It should be noted that this program assumes the earth to be spherical 
rather than oblate as it actually is. However, by assuming a spherical 
earth, one introduces maximum error distances of about 11.5 nautical 
miles. For many purposes, this error can be ignored. 



II. MATHEMATICAL ANALYSIS 

Figure 1 shows the intersection geometry. The satellite is located 
at point SAT in a circular, near-equatorial orbit which is a distance 
DIST from center of the earth. Point P is the point of intersection of 
satellite antenna boresight (or beam center) and the earth. AL is the 
arc from the point on earth directly below the satellite (the subsatellite 
point) to P and is given by: 

AL = Cos" 1 [Cos ( LONCTR) Cos ( LACTR + DELT)] (1) 

Where LATCTR is latitude of P, LONCTR is the longitude of P relative to 
the subsatellite point, and DELT is the instantaneous declination angle 

formed by a vector from satellite to the center of the earth and the 

22,3t>0 

equatorial plane. At DIST = 30,330 statute miles (geostationary equa- 
torial orbit when DELT =0), the maximum arc, AL, permissible for the 
point to be seen by the satellite is 81.3°. When AL is computed to be 
larger than 81.3°, the point is over the horizon seen by the satellite 
and any attempt at finding the locus of intersection would produce 
meaningless results. Therefore, the computation stops at this point. 

The next step is to calculate the vector, RS, extending from the 
satellite to point P. A right-handed coordinate system is constructed 
with the origin at SAT, the positive Y axis extending through the earth 
center and the Z axis in the plane defined by SAT and the north and south 
poles. The vector RS is given in terms of coordinates along these axes: 

RS = RE cos (LATCTR) si n( LONCTR) i_ + 

[DIST cos (DELT) - RE cos(LATCTR) cos(LONCTR)] j_ + 

[DIST sin(DELT) + RE sin (LATCTR)] k (2) 

Where RE is radius of earth, and i_, j_ and are unit vectors in 
the positive x, y, and z directions, respectively. 

Pitch and roll angles are then defined as in Figure 2. These are 
given by: 







PITCH = Sin' 1 


(3) 


RS. 


/ — 2 T 

VRS i + RS, 


ROLL 



(4) 


A new coordinate system is defined by rotating the old system about 

I 

its z-axis by an angle equal to PITCH, and then about the resultant x-axis 
by an angle equal to ROLL. 

A vector L[ is then defined as lying in, a plane perpendicular to 
vector RS^. The origin of vector is at the point of intersection of 
RS and that plane, see Figure 3. 

The equation defining U_ is: 

U = [cos (BETA) i* + sin(BETA) k'] [A“ 2 cos 2 (BETA) + 

B- 2 sin 2 (BETA)] -Js (5) 

i 

where i_' and Ik* are unit vectors in the new coordinate system. This 
equation is a parametric equation of an ellipse with BETA as the parameter 
and A and B corresponding to the semi -major and semi -mi nor axes, 
respectively. 

A ray on the surface of the beam can now be generated by a vectorial 
addition of RS and U_. By incrementing BETA from 0° to 360° the entire 
outer surface of the beam can be generated. 

Let a vector Mn be a surface generator vector. The equation defining 
Mn in the new coordinate system is: 

Mn = [ | RSj cos (ROLL) sin(PITCH) + |U| cos(BETA) cos(PITCH) 

- |U| sin(BETA) sin(ROLL) sin(PITCH)] 1_' 

+ [ | RSj cos (ROLL) cos (PITCH) + |U| cos(BETA) sin(PITCH) 

- |U| sin(BETA) sin(ROLL) cos(PITCH)] y 

+ [ | RS | si n( ROLL) + |U| sin(BETA) cos ( ROLL)] k' (6) 
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In the old coordinate system: 


Mn = [Mn. ,] i + [Mn., cos(DELT) + Mn. , sin(DELT)] £ 

I J K 

+ [-Mn - , sin(DELT) + Mn. , cos(DELT)] k (7) 

Pitch and roll angles are defined for each Mn as follows: 

Mn i , 

PITCH N = Sin" 1 (8) 

/ 2 2 ~ 

/v/Mn . , + Mn^. , 

M V 

ROLL N = Tan” 1 — ■ (9) 


The earth radius vector to the point of intersection of Mn is: 

) 

RE = [Mn cos(ROLLN) sin(PITCHN)] I 

+ [DIST cos(DELT) - Mn cos(ROLLN) cos(PITCHN)] J 
+ [Mn sin(ROLLN) - DIST sin(DELT)] K (10) 


Where Mn in the length of Mn and and J< are unit vectors in an 

earth centered coordinate system whose positive z axis extends through 
the north pole and positive x axis intersects the 0° meridian. 

The latitude and relative longitude coordinates of the point of 
intersection of Mn are found by: 


LAT = Tan 


RE. 


-l 


/ 2 2 

a7RE t + RE, 


( 11 ) 


LON = Sin' 


RE 


I 


^/REj + REj 


( 12 ) 
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The actual longitude is found by adding the above longitude to the 
longitude of the subsatellite point. 

The relative beamwidth for any signal level is found using a beam- 
width conversion chart. The particular chart used in this program is taken 
from Microwave Engineers Handbook^, (see Figure 4). The actual con- 
version in the program is done by using linear interpolation between the 
appropriate two consecutive points from the following set: 

(.1, .18) (.2, .26) (.5, .4) (1., .56) (1.5, .7) (3., 1.) (5., 1.27) (10., 1.7) 


These points were chosen such that the graph segment between any two 
consecutive points is approximately linear. 

The above discussion assumes that the entire antenna beam intersects 
the earth. The case when the boresight location is near enough to the 
horizon that a portion of the beam passes the earth is considered next. 

Define an angle, Bn formed by the vector Mn and the vector from the 
satellite to earth center. The angle at which Mn is tangent to the earth 
is given, from the law of sines, as: 


B 


max 


Sin" 1 


RE 

DIST 


(13) 


When the angle B„ is larger than B the earth radius vector RE 

must be computed to the point of tangency, i.e. the horizon seen by the 
satellite. The vector RE^ is then defined by: 

RE = RE [cos(B max ) sin(TAUN)] I 

- RE [cos(DELT) sin(B max ) 

+ sin(DELT) cos(B max ) cos(TAUN)] J 

- RE [sin(DELT) sin(B max ) 

- cos(DELT) cos(B m 1 cos(TAUN)] K (14) 

When TAUN is the tilt angle defined by: 


Mn k 

TAUN = Cos' 1 = 


V Mn i + Mn k 


(15) 
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and RE is the radius of the earth. The latitude and relative longitude 
coordinate are found as in equations (11) and (12). 


III. USAGE 

The inputs to the program are: 

LONSS - longitude of point on earth directly below satellite (°E) 

LONCT - longitude of bores ight intersection (°E) 

LATCT - latitude of boresight intersection (°N) 

INCR - increments of angular parameter BETA, i.e. number of points 
plotted = 3607INCR 

DELT - instantaneous declination of satellite - earth center vector from 
plane of equator (°below equatorial plane) 

THETA - orientation angle of elliptical beam about beam center (measured 
positive counter-clockwise) 

BW1 - beam vertex angle in plane defined by beam vertex and major axis 
of beam cross section at point 3 db below beam center 
BW2 - beam vertex angle in plane defined by beam vertex and minor axis 
of beam cross section at point 3 db below beam center 
L - number of signal levels to be plotted 

DB(I) - signal levels to be plotted (measured in db below beam center) 

It should be noted that all inputs except L and DB ( I ) are in degrees 
rather than radians. Longitude inputs can be expressed as degrees east 
of 0° or negative degrees west of 0°. The output will agree with the 
input. 

Each input data set will consist of three cards. The format for 
the data is: 

Card #1 

column 
1 - 9 

10 - 18 
19 - 27 
28-36 
37 - 45 


LONSS 

LONCT 

LATCT 

INCR 

DELT 
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Card #2 


1 

10 

19 

28 


9 

18 

27 

29 


THETA 

BW1 

BW2 

L 


Card #3 


1 

8 

15 


7 

14 

21 


DB(I) 


Every input, with the exception of L, must contain a decimal point. 
The input L must not contain a decimal point and must be right justified 
in columns 28-29, i.e. single digit input must be in column 29. The 
inputs DB(I) must be input in increasing order with the signal level 
closest to beam center listed first and the level farthest from beam 
center listed last. 

The output of the program consists of two parts, the printout and 
the plot. The printout contains, for each set of input data, a list 
of the following parameters: 

SUB SAT LONG - longitude of point directly below satellite 
BOR SGHT LONG - longitude of bores ight intersection 
BOR SGHT LAT - latitude of boresight intersection 

DECLINATION - instantaneous declination of satellite-earth center vector 
from equatorial plane 

MIN BMWDTH - beamwidth along minor axis of beam cross section at -3 db 
(half-power) level 

MAX BMWDTH - beamwidth along major axis of beam cross section at -3 db 
(half-power) level 

ORIENTATION - orientation angle of beam about beam center 

ELEVATION - angle formed by vector RS and a plane tangent to earth at 

boresight intersection 

The remainder of the printout gives the maximum and minimum beamwidth 
and a listing of the coordinates of the locus of intersection for each 
signal level. 
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Each set of input data is given a data set number. This number 
appears on the printout and the plot for each data set. This facilitates 
matching of plots with printout when more than one data set is run. 

A modification of this program has recently been developed to plot 
the coverage of a multi -beam satellite. This new program computes the 
locus of intersection of a number of sets of input data as does the original. 
However, the modification plots all of the intersection loci on one set of 
axes so that uncovered areas and overlapped areas are immediately obvious. 

The major changes to the original program are: 

1. The plotter tape is opened and the axes are drawn before the 
main calculation begins. 

2. The axes are scaled once at the beginning of the program rather 
than being scaled for each data set. 

3. The computed coordinates are checked to see that they do not 
extend beyond the limits of the axes. 

4. The plot origin is not reset for each data set. 


IV. TYPICAL EXAMPLES 

Figures 5-9 present the area coverage plots for input values shown 
in Table 1-. For each case, 3, 5 and 10 dB level contours are plotted. 

The inputs represent a wide variety of cases--satellites positioned in 
circular, equatorial and geosynchronous orbits; satellites positioned in 
slightly inclined stationary orbits; and area coverage at small inclination 
angles. 

Figure 1 0 shows the coverage provided by a 4-beam satellite positioned 
at 120°W longitude. All plots use the same set of axes. The input values 
for various beams are given in Table 2. All four beams are designed to 
provide complete coverage to the United States— beam 1 covers Hawaii, 
beam 2 covers Alaska, and beams 3 and 4 provide coverage to the other 
48 states. 


0 
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PROGRAM LISTING 
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